Skip to contents

Lake Como, nestled in the Lombardy region, is recognized as a crucial freshwater ecosystem esteemed for its biodiversity and ecological significance. Suppose that environmental studies have suggested the presence and spreading of a concerning pollutant within the lake’s waters. The emergence of this pollutant could raise significant concerns about its potential threat to the ecological balance of the aquatic environment and necessitates considerations about its long-term environmental stability. This hypothetical presence prompts a thorough examination to understand its nature, source, and potential implications on Lake Como’s delicate ecosystem.

Assume that several monitoring stations have been strategically established around Lake Como to track and measure pollutant concentrations in its waters. These stations play a pivotal role in continuous surveillance, offering crucial data to assess pollutant levels and their potential impacts on the lake’s ecosystem.

The following interactive figure shows the comprehensive map of Lake Como and includes georeferenced black points denoting the strategic placement of monitoring stations across its expanse. The red points denote potential ‘hotspots’ where, in the hypothetical scenario under investigation, abnormal concentrations of the pollutant have been recorded.

In modeling the spread of the pollutant, it’s crucial to consider the presence of a tributary in the northern region of the lake. Consequently, a transport term has been identified to represent the movement of the pollutant carried by the inflow to other areas within Lake Como. Furthermore, non-homogeneous Dirichlet boundary conditions are applied to boundary regions near the stations that recorded anomalies in the pollutant concentration. Conversely, homogeneous Dirichlet boundary conditions are applied to boundary regions unaffected by the anomalies.

Steady-state solution

In this section, we focus on solving the steady-state problem. Let \(\Omega\) represent our domain of interest. We assume that the concentration of the pollutant, denoted by \(u\), spreads across \(\Omega\) according to the following partial differential equation (PDE):

\[ \begin{cases} -\mu \Delta u + \boldsymbol{\beta} \cdot \nabla u = 0 \qquad & in \ \Omega \\ u = g|_{\partial \Omega} \qquad & on \ \partial \Omega, \end{cases} \] where \(\mu \in \mathbb{R}\), \(\boldsymbol{\beta} \in \mathbb{R}^2\) are the diffusion coefficient and the transport coefficient, respectively. The function \(g\) used to enforce the boundary conditions consists of the sum of two Gaussian-like functions centered at the locations of the stations that detected anomalous concentrations of the pollutant.

The upcoming sections will illustrate how to obtain a discrete solution of the problem using femR package.

library(sf, quietly = T)
library(mapview, quietly = T)
library(femR, quietly = T)
library(rmapshaper, quietly = T)

lake_bd <- read_sf("data/deims_sites_boundariesPolygon.shp")
lake_bd_simp <- rmapshaper::ms_simplify(lake_bd, keep= 0.5, 
                                        keep_shapes=T)
st_crs(lake_bd_simp) <- 4326

The femR package includes utilities to read a sfc (Simple Features Collection) object, such as the boundary of the Lake, and generating meshes by specifying the maximum_area of each triangle and the minimum_angle between adjacent sides of the triangles.

# create the physical domain
domain <- Domain(st_geometry(lake_bd_simp)) 
# create the mesh
mesh <- build_mesh(domain, maximum_area = 0.25e-5, minimum_angle = 20)

Furthermore, the Mesh class overloads the st_as_sfc method, enabling visualization of the mesh object on an interactive map using the mapview package, as demonstrated in the following code snippet.

mesh_sf <- st_as_sfc(mesh)

map.type <- "Esri.WorldTopoMap"
mapview(lake_bd_simp, col.regions = "white", lwd = 2, 
          alpha.regions=0.2, legend=F, map.type=map.type) +
mapview(mesh_sf, col.region="white", alpha.regions=0.2, legend=F,  alpha=0.4,
        lwd=0.4, map.type=map.type) 

The next step involves defining the FunctionSpace where to find the solution of the PDE. By default, first-order finite elements are considered during the construction of Vh. Then, the solution of the problem is instantiated as an element of Vh.

Vh <- FunctionSpace(mesh)
u <- Function(Vh)

Then, we define the differential operator in its strong formulation as follows:

mu <- 0.001
beta <- c(-4.091349, -4.493403)/20
Lu <- -mu*laplace(u) + dot(beta, grad(u))

Both the forcing term of the differential problem and the boundary condition can be defined either as standard R functions or as simple numerical values, particularly when representing constant values.

# Dirichlet boundary conditions 
g <- function(points){
  res <- matrix(0, nrow=nrow(points), ncol=1)
  pts_list <- lapply(as.list(as.data.frame(t(points))), st_point)
  dists2 <- st_distance(st_geometry(sources),
                        st_sfc(pts_list, crs=st_crs(sources))) 
  dists2 <- matrix(as.numeric(dists2), 
                   nrow=length(st_geometry(sources)), ncol=nrow(points))/(2.5*10^3)
  for(i in 1:length(sources))
    res <- res + Q*exp(-dists2[i,]) 
  return(res)
}

Finally, we build a Pde object passing the differential operator L, as first parameter, and the forcing term, which is a constant in this specific case, as second parameter:

# build PDE object
pde <- Pde(Lu, 0.)
# set boundary conditions
pde$set_boundary_condition(fun= g, type= "dirichlet")

We can compute the discrete solution of the problem calling the solve method:

pde$solve()

Finally, it is possible to visualize the computed solution on an interactive map using the mapview package as the following code snippet shows.

coeff <- apply(mesh$get_elements(), MARGIN=1, 
               FUN= function(edge){
                     mean(u$coeff[edge])
         })
U <- st_sf(data.frame(coeff = coeff), geometry= mesh_sf)
mapview(U, alpha.regions=0.8, alpha = 0, map.type=map.type)

Time-dependent solution

Now, we aim to estimate the evolution of pollutant concentrations in Lake Como. We consider a Gaussian-like function \(g\) from the previous section, appropriately modified to account for the presence of time. This function will be used to define both the initial condition for the pollutant concentration and to impose the boundary conditions. Therefore, the partial differential equation modeling the evolution of the pollutant concentration in Lake Como is as follows:

\[ \begin{cases} \frac{\partial u}{\partial t} -\mu \Delta u + \boldsymbol{\beta} \cdot \nabla u = 0 \qquad & in \ \Omega \times (0,T) \\ u = g & in \ \Omega, \ t=0 \\ u = g|_{\partial \Omega} \qquad & on \ \partial \Omega \times (0,T), \end{cases} \]

The following window wraps all the steps needed to estimate the evolution of the pollutant concentration relying on femR package.

time_interval <- c(0., 0.3)
times <- seq(time_interval[1], time_interval[2], length.out=50)
# function space
Vh <- FunctionSpace(mesh %X% times)
# PDE solution
u <- Function(Vh)
# differential operator
Lu <- dt(u) - mu*laplace(u) + dot(beta, grad(u))
# dirichlet bc
g <- function(points, times){
  res <- matrix(0, nrow=nrow(points), ncol=length(times))
  pts_list <- lapply(as.list(as.data.frame(t(points))), st_point)
  dists2 <- st_distance(st_geometry(sources),
                        st_sfc(pts_list, crs=st_crs(sources))) 
  dists2 <- matrix(as.numeric(dists2), 
                   nrow=length(st_geometry(sources)), ncol=nrow(points))/(2.5*10^3)
  for(i in 1:length(sources)){
      res[,1:length(times)] <- res[,1:length(times)] + Q*exp(-dists2[i,])
  }
  return(res)
}
# initial condition
u0 <- g(mesh$get_nodes(), times[1])
# build PDE object
pde <- Pde(Lu, 0.)
# set initial conditions
pde$set_initial_condition(u0)
# set boundary conditions
pde$set_boundary_condition(fun= g, type= "dirichlet")
# compute the discrete solution
pde$solve()
# plot 
plot(u)